  ' Rice pdf enter gamma, ga and sigma,sb   from csc pdf fit 7/6/86
  
  'Q basic version, including graphics for Mac
  
  DIM e(500),WR(4,500),I(12),j(12)
  
  'input is e(m) and the output is WR(j,m) the rician pdf for
  'of gamma
  
  PI = 4 *  ATN (1)
  mmax = 500
  m1 = 1
  a = .004
  fw = 1
  
  'following values of gamma give a set for graphs
  
  gam(0) = 0: gam(1) = 1: gam(2) = 2 :gam(3) = 5: gam(4) = 10
  
  'CALC MOD BESSEL FUNCTION I0(X).USE ABRM AND STEGUN POLYNOM.
  'the coefficients follow
  
 I(0) = 3.51562:I(1) = 3.089942:I(2) = 1.206749
 I(3) = .265973:I(4) = .0360768:I(5) = .0045813
  j(0) = .39894228#
  j(1) = .01328592#:j(2) = .00225319# :j(3) = -.00157565#
  j(4) = .00916281#:j(5) = -.02057706#
  j(6) = .02635537#:j(7) = -.01647633#: j(8) =.00392377#
 
 PRINT" Rician pdf" 
  PRINT "input sigma =";:INPUT sb
  
10 ' start of rice pdf calc
 'the parameters are ga and sb
 'Gamma goes to the sub as ga
 'Sigma is sb
 
 'Rice is computed as a function of e(m), proportional to x in the sub.
 
   FOR j = 0 TO 4
       PRINT "gamma =";gam(j)
       FOR m = m1 TO mmax
        e(m) = a*m
        ga = gam(j)
      GOSUB 1000
       NEXT m
   NEXT j
   
   REM  	SCREEN DIMENSIONS
     XL = 460
     YL = 260
 
REM   	SET SCALES

 X0 = 20
     XS = (XL - X0) / e(mmax): REM  X(NM) IS MAXIMUM VALUE OF X
     Y0 = 270: REM  TO PUT Y=0 IN MIDDLE
     YS = 240/2: REM  THIS SETS THE AMPLITUDE FACTOR.
 
REM	TOOL BOX CALLS REQUIRE INTEGERS. % INDICATES INTEGER  
REM  	CALCULATE X% AND Y% AND THEN PLOT TO X1% AND Y1%.  

 CLS         : REM CLS clears the screen
 PICTURE ON  : REM PICTURE ON puts screen graphics in storage.
 SHOWPEN     : REM SHOWPEN also puts graphics on the screen 
FOR j = 0 TO 4
   FOR n = 1 TO mmax - 1
          x% =  INT (XS * e(n) + X0)
     X1% =  INT (XS * e(n+1) + X0)
          y% =  INT (Y0 - YS * WR(j,n))     : REM plot + up.
          Y1% =  INT (Y0 - YS * WR(j,n + 1))
          LINE (x%,y%)-(X1%,Y1%)         : REM draw line
   NEXT n
NEXT j

 x% = INT (X0) : X1% = INT (XS*e(mmax) + X0)
 y% = Y0
 LINE (x%,y%) - (X1%,y%)            :REM draw axis
 LINE (x%,Y0) - (x%,Y0-2*YS)
 
 FOR n = 1 TO 20
  x% = X0
 y% = INT (Y0-YS*n/10)
 LINE (x%,y%) - (x%+2,y%)
NEXT n
 
REM	PUT TICS ON THE X-AXIS

 FOR n = 0 TO mmax STEP mmax/10
      x% =  INT (XS * e(n) + X0)
     y% = INT (Y0 )                 :REM make tics
      LINE (x%,y%) - (x%,y%-2)            :REM draw tics
     CALL MOVETO (x%-10,y%+20) : PRINT e(n);
 NEXT n
 
 'CALL MOVETO (20,290) : PRINT "sb=";sb
     
 PICTURE OFF : REM PICTURE OFF ends graphics operations.
INPUT Q$
IF Q$ = "ng" GOTO 100
 CALL MOVETO (20,280)   
   PRINT " input 'y' to make a file "
    INPUT Q$
 IF Q$ <> "y" GOTO 100

 pic$ = PICTURE$          :REM PICTURE$ is name of stored picture.
 
 CALL MOVETO (50, 25)     :REM name the file 
 PRINT  "I've got the picture in pic$ ("; LEN (pic$); ")"
 pictFile$ = FILES$ (0, "Enter name for PICT file:")
 PRINT "PICT file name is:"; pictFile$
 
REM	SAVE FILE IN 'PICT' FORMATE.

 OPEN pictFile$ FOR OUTPUT AS #1
 
REM	FOR-NEXT LOOP MAKES A HEADER FOR PICT FILE FORMATE.

 FOR I = 1 TO 512 : PRINT  #1, CHR$ (0); : NEXT
 
 PRINT  #1, pic$
 
 CLOSE  :REM the picture 'pic$' is stored as a text file.
 
REM	CHANGE THE FILE TYPE FROM TEXT TO PICT

 NAME pictFile$ AS pictFile$, "PICT"
 
REM	USE MacDraw TO READ THE FILE. THEN, 
REM	IT CAN BE SAVED AS A MacDraw DRAWING.
 
100  CLS    :REM clear screen and end


   
   INPUT Q$
   
   END
   
1000 'rice pdf subroutine

  XE = ((1 + ga) * e(m)^2 + ga * sb) / sb
  x = 2 * e(m) *  SQR (ga * (1 + ga)) /  SQR (sb)
  XC = 2 * e(m) * (1 + ga) / sb
   
   IF x < 3.75 THEN 

 T = x / 3.75

 T2 = I(0) * T ^ 2
 T3 = I(1) * T ^ 4
 T4 = I(2) * T ^ 6
 T5 = I(3) * T ^ 8
 T6 = I(4) * T ^ 10
 T7 = I(5) * T ^ 12
 I0 = 1 + T2 + T3 + T4 + T5 + T6 + T7
      WR(j,m) = fw * XC *  EXP ( - XE) * I0
   ELSE
 T = x / 3.75
 
 T1 = j(1) * T ^  -1
 T2 = j(2) * T ^  -2
 T3 = j(3) * T ^  -3
 T4 = j(4) * T ^  -4
 T5 = j(5) * T ^  -5
 T6 = j(6) * T ^  -6
 T7 = j(7) * T ^  -7
 T8 = j(8) * T ^  -8
 I1 = j(0) + T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8
      WR(j,m) = fw *XC *EXP(x - XE)*I1/ SQR(x)
   END IF
   
   IF YM < WR(j,m) THEN YM = WR(j,m)
   
   RETURN
   
  

